//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs

rename cum_avg_daycare_36m_sum c
rename parenting_ages13_std    p

// childcare
// average for singleton control
summ c if twin == 0 & tg == 0
local c_s0_t0 = r(mean)
di round(r(mean),.01)

// average for twin pair 
summ c if twin == 1 & tg == 0
di round(r(mean),.01)
local c_s1_t0 = r(mean)
di 2*round(r(mean),.01)

// parenting
// average for singleton control
summ p if twin == 0 & tg == 0
di round(r(mean),.01)
local p_s0_t0 = r(mean)

// average for twin pair 
summ p if twin == 1 & tg == 0
local p_s1_t0 = r(mean)
di round(r(mean),.01)
di 2*round(r(mean),.01)

// ratios
local  pc_ratio_s0_t0 =  `p_s0_t0' /  `c_s0_t0'
di round(`pc_ratio_s0_t0',.01) 
matrix pc_ratio_s0_t0 = `pc_ratio_s0_t0'

local pc_ratio_s1_t0 =  `p_s1_t0' /  `c_s1_t0'
di round(`pc_ratio_s1_t0',.01) 
matrix pc_ratio_s1_t0 = `pc_ratio_s1_t0'

// difference in ratios
di round(`pc_ratio_s1_t0' - `pc_ratio_s0_t0',.001)
matrix pc_dratio_t0 = `pc_ratio_s1_t0' - `pc_ratio_s0_t0'

// difference in ratios, 2 
di round(`pc_ratio_s1_t0' - `pc_ratio_s0_t0',.001)
matrix pc_ddratio_t0 = `pc_ratio_s1_t0' - 2*`pc_ratio_s0_t0'

// difference in control parenting
di round(`p_s1_t0' - `p_s0_t0',.001)
matrix pdiff_t0 = `p_s1_t0' - `p_s0_t0'

// diff double 
di round(`p_s1_t0' - `p_s0_t0',.001)
matrix pdiff_t0 = `p_s1_t0' - `p_s0_t0'

foreach var in pc_dratio_t0 pdiff_t0 { 
	matrix `var'_B = [.]
}

// bootstrap
foreach b of numlist 1(1)$bootstraps {
	preserve
	bsample, strata(bwg site tg)
	// childcare
	// average for singleton control
	summ c if twin == 0 & tg == 0
	local c_s0_t0 = r(mean)

	// average for twin pair 
	summ c if twin == 1 & tg == 0
	local c_s1_t0 = r(mean)

	// parenting
	// average for singleton control
	summ p if twin == 0 & tg == 0
	di round(r(mean),.01)

	// average for twin pair 
	summ p if twin == 1 & tg == 0
	local p_s1_t0 = r(mean)
	
	// ratios
	local  pc_ratio_s0_t0 =  `p_s0_t0' /  `c_s0_t0'
	local  pc_ratio_s1_t0 =  `p_s1_t0' /  `c_s1_t0'

	// difference in ratios
	matrix pc_dratio_t0_`b' = `pc_ratio_s1_t0' - `pc_ratio_s0_t0'
	matrix pc_dratio_t0_B = [pc_dratio_t0_B \ pc_dratio_t0_`b']
	
	// difference in control parenting
	matrix pdiff_t0_`b' = `pc_ratio_s1_t0' - `pc_ratio_s0_t0'
	matrix pdiff_t0_B = [pdiff_t0_B \ pdiff_t0_`b']
	restore
}

// inference
foreach var in pc_dratio_t0 pdiff_t0 { 
	matrix `var'_B = `var'_B[2...,1...]
	matrix colnames `var'_B = `var'_B
	
	svmat   `var'_B, names(col)
	summ    `var'_B 
}

// tests hypotheses 
// difference of ratios 
replace pc_dratio_t0_B  = pc_dratio_t0_B - pc_dratio_t0[1,1]
gen     pc_dratio_t0    = pc_dratio_t0[1,1]

gen     indpc_dratio_t0 = 0 if pc_dratio_t0_B !=. 
replace indpc_dratio_t0 = 1 if pc_dratio_t0_B > pc_dratio_t0 & pc_dratio_t0 !=. & pc_dratio_t0_B !=. 
summ    indpc_dratio_t0

// difference in parenting
replace pdiff_t0_B  = pdiff_t0_B - pdiff_t0[1,1]
gen     pdiff_t0    = pdiff_t0[1,1]

gen     indpdiff_t0 = 0 if pdiff_t0_B !=. 
replace indpdiff_t0 = 1 if pdiff_t0_B < pdiff_t0 & pdiff_t0 !=. & pdiff_t0_B !=. 
summ    indpdiff_t0
